Covariant Kinetic Theory with an Application to the Coma 

Cluster 

Brandon Wolfe 1 and Fulvio Melia 1,2 

1 Physics Department, The University of Arizona, Tucson, AZ 85721 

2 Steward Observatory, The University of Arizona, Tucson, AZ 85721 

ABSTRACT 

In this paper, we introduce a novel solution to the covariant Landau equation 
for a pure electron plasma. The method conserves energy and particle number, 
and reduces smoothly to the Rosenbluth potentials of non-relativistic theory. In 
addition, we find that a fully relativistic plasma equilibrates in only 1/ 100th 
of a Spitzer time — much faster than in the non-relativistic limit — a factor of 
significant import to situations in which distortions to a Maxwellian distribu- 
tion are produced by anomalous methods of acceleration. To demonstrate the 
power of our solution in dealing with hot, astrophysical plasmas, we use this tech- 
nique to show that one of the currently considered models — continuous stochastic 
acceleration — for the hard X-ray emission in the Coma cluster actually cannot 
work because the energy gained by the particles is distributed to the whole plasma 
on a time scale much shorter than that of the acceleration process itself. 

Subject headings: acceleration of particles — galaxies: clusters: individual (Coma) 
- plasmas — radiation mechanisms: non-thermal — relativity — X-rays: galax- 
ies 



1. Introduction 

The need for a self-consistent, practical transport theory that can handle dynamic pop- 
ulations of relativistic particles extends across many disciplines in physics. 

Most obviously, astrophysical plasmas in the presence of strong gravitational and/or 
magnetic fields are often out of thermal equilibrium when the relevant dynamic time scale 
(e.g., associated with the process of accretion onto a compact object) is short compared to the 
time required for the particles to interact internally with each other. Situations in which this 
may occur include solar flares (Petrosian and Liu 2003), accretion onto supermassive black 
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holes in the nuclei of active galaxies (see, e.g., Melia and Falcke 2001), and re-acceleration 
of already relativistic particles in shocks produced during the merger of galaxies (Brunetti 
et al. 2004). In all cases, the resultant photon spectrum produced by the energized particles 
is most simply explained in terms of a non-Maxwellian distribution. 

Less directly, an argument can be made that the two dominant particle species (say, 
electrons and protons in a fully ionized Hydrogen plasma) separate in energy space, leading 
to a situation in which the ion temperature is larger than that of the electrons (see, e.g., 
Shapiro, Lightman, and Eardley 1976). However, whether such a two-temperature plasma 
may be maintained against collisional re-equilibration is still actively debated, and may be 
critical to the question of how a black hole accretes from its environment. For example, there 
is now very good evidence that the gas orbiting the supermassive black hole, Sgr A*, at the 
center of our Galaxy attains a temperature in excess of 10 10 K (Liu and Melia 2001), but it 
is still not understood whether this optically thin plasma maintains a single temperature. In 
addition, nonthermal emission in Sgr A* (see, e.g., Ozel, Psaltis, and Narayan 2000; Yuan, 
Quataert, and Narayan 2003) has been invoked to explain its low-frequency radio spectrum 
in accretion-dominated models. 

These clear astrophysical applications should be aided by recent experiments with laser- 
heated plasmas, which are now probing the energy regimes of interest. Phenomena unique to 
relativistic distributions, such as their tendency to limit heat flux (Bell et al. 1981; Rickard 
et al. 1989), or the ability of rapid acceleration to create bi-Maxwellian electron distribu- 
tions (Guethlein et al. 1996) have now been observed in laboratory plasmas. Experimental 
constraints arising from both the astrophysical and laboratory contexts, acting in concert, 
should yield a solid, well-tested theory. 

Astrophysicists have had a theory for relativistic energy transport for some time. Dewar 
et al. (1979) were followed by Dermer and Liang (1989) in calculating the energy exchanged 
between particles during a relativistic binary collision. Nayakshin and Melia (1998) used 
these coefficients as the basis of a Fokker-Planck equation in energy. But, as Dewar et al. 
note, "virtually all of the quantities that were invariant in the non-relativistic treatment 
cease to be so in the relativistic formulation. Much care (must) be taken to show how the 
various physical quantities... transform." The resulting theory is not as transparent as one 
would like. 

In contrast, the strategy we adopt here follows Beliaev and Budker (1955) in formulating 
kinetic theory in a covariant form valid in any frame of reference. This approach is simpler, 
since the frame invariants are obvious. It is also general: we present a solution that handles 
anisotropic or radically out-of-equilibrium distributions, and can be used to find electric and 
thermal conduction coefficients as well as energy exchange. This approach also conveniently 
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and elegantly reduces to the well-known Rosenbluth (1957) theory in the non-relativistic 
limit. 

Relativistic distributions obey qualitatively different physics, since the covariant phase 
volume must incorporate energy to form an eight-dimensional phase space. But the covariant 
theory is also quantitatively sufficiently different to shape observations. We show here that 
a fully relativistic distribution equilibrates in a tenth to a hundredth the time given by 
non-relativistic theory. The heat flux in a relativistic distribution, relative to that for the 
non-relativistic case, is inhibited by a factor of 1CT 2 — 1CT 1 . This is due to a drift velocity 
which approaches the speed of light asymptotically. Thus, situations where a nonthermal 
tail is likely to be present (see, e.g., Blasi 2000; we will examine this situation in greater 
detail below) require a covariant formulation. 

Solutions to covariant kinetic theory have existed for some time, too. Braams and 
Karney (1987) made a spherical expansion of the kinetic kernel and found relativistic con- 
ductivities (1989), but only in very limited circumstances. This was followed by attempts at 
solution via a Chapman-Enskog expansion (Mohanty and Baral 1996), Legendre expansion 
(Honda 2003), Grad expansion (Muronga 2004), and numerical integration (Shoucri and 
Shakarofsky 1994). 

But, to our knowledge, none of these solutions is six dimensional — the ability to solve 
three dimensional integrals like (Eq. 37 below) has only existed for a short time (Hahn 
2004). Manheimer et al. (1997) have called the self-consistent evaluation of even the simpler 
non-relativistic coefficients "unfeasible". 

Many approach the Fokker-Planck equation with implicit differencing, which is relatively 
difficult to modify with additional forces (e.g., in a particle-in-cell simulation). Instead, we 
follow the original suggestion of Jones et al. (1996) and the subsequent application of Habib 
et al. (2000) to model collisions as Brownian motion. A single particle's motion is observed at 
discrete times, and the affect of multiple collisions accumulated over this time is represented 
as dynamical friction and noise. The distribution is approximated by some small fraction 
of the actual number of particles, all moving according to self-consistent kinetic forces. The 
simulation thus solves single-body statistical dynamics, many times. 

The desire for a covariant theory that is simple to apply and reduces naturally to the 
non-relativistic version is far from academic. As a concrete — and important — example, we 
show that our theory eliminates the possibility that stochastic acceleration is responsible 
for the hard X-ray emission in the Coma cluster. Before that, however, we present the 
relativistic kinetic theory in § 2, and then describe the method of solution in covariant form 
in § 3. The application will be made in § 4. 
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2. Relativistic Kinetic Theory: Background 

We begin by discussing the development of the relativistic kinetic equation somewhat 
pedagogically, given the fact that its usage in the analysis of high-energy plasmas in astro- 
physics has been infrequent at best. Landau (1936) first approached the problem by showing 
that the Boltzmann collision integral can be written as a flux in velocity space, 

thereby writing the Boltzmann equation as a conservation equation for the single particle 
phase space distribution, 

^/ + {ff,/} = V-S, (2) 

where / is the single particle phase space number density. The Poisson bracket {H, /} 
represents a small change in the shape of a volume of phase space (which nevertheless 
conserves number), the term y . S a flux of particles out of it. But what is S7 To begin 
with, note that the number of collisions occurring in unit time between a particle with 
velocity v and particles with velocity v' in the range d 3 v ' is 

™/(v)/(v') d 3 q d 3 v> , (3) 

where q is the amount of velocity exchanged. Here, w, the probability that v will be taken 
into v + q, can be expressed evenly as w(y + q/2, v' — q/2;q). This ensures that w is 
symmetrical with respect to q — i.e., that detailed balance is maintained. 

The particle flux F = F + — F~ is found by subtracting the number of particles entering 
a volume of velocity space from the left, 

F~ = I d 3 q j d s v' r wf(v)f'(v') dv a (4) 

from that exiting to the right, 

F + = [ d 3 q f d 3 v' f a wf(v + cO/V - q) dv a . (5) 

This leads to an argument /(v)/'(v') — /(v + q)/'( v ' — q) which can be expanded about 
small q, giving the required S, the so-called Landau collision integral, 



c = —s 9 V 



/(# - /'(V) s/W 



dv' p dvp 



B aP d 3 v , (6) 
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where w is expressed in terms of the collision cross section, 

w d 3 q — |v — v'| da , 

and where 



B a( s = \ J Q»q^ - v'| da 



(7) 
(8) 



Finally, since momentum conservation dictates that the change q is perpendicular to the 
relative velocity |v — v'|, it must be that 

B af3 (vp -v' (3 ) = 0. (9) 

And since B a/3 can only depend on the vector |v — v'|, the tensor must take the form 



B a/ 3 = -B 



) a 



(v - v') 2 



where 



B = B aa = i J q 2 \w - v'| da . 
Thus, integrating Equation (11) with a Rutherford cross-section for da, one gets 



B af 3 — 

where U(v, v') is the collision kernel, 

U(v,v') = 



e 2 e' 2 



lnAU(v,v') , 



(10) 



;n) 



(12) 



v — v' 




- (v a - 


-v' a )(v{, -v' p ) 




v — v' 


CO 



(13) 



and a, (3 are Cartesian indices. In what follows, we define A = ne 2 e' 2 /8irelm 2 In A and 
I = a/ mjlkT for succinctness. In c.g.s. units, A = 87re 2 e /2 nlnA/m 2 . 

With one integration by parts, the collision integral can be written 



C 



_d_ 

0v n 



D 



0/(v) 



- F a /(v) 



where 



and 



D a p = Aj U(v,v')/(v') rfV, 



F a = Aj2 



d 



dv 



-U(v,V) 



(14) 

(15) 
(16) 
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which is the Fokker-Planck equation seen commonly in the west. Equations (15) and (16) 
can be expressed in a simpler form by observing that the kernel U(v, v') obeys the relations 



and 



d 

dv' a 



U n = 



<9 2 |v — v'| 

dv a dv[3 

-2d\v- V]- 1 



(17) 
(18) 



thus allowing the coefficients to be expressed as 

d 2 



D aP = —A- 



and 



Fa=-A„ 

on 



dVadV/3 

d 



f\v-v'\d 3 v' 



2 J /|v- v'l -1 d 3 v' 



(19) 
(20) 



The advantage of this phrasing (Rosenbluth et al. 1957) is that the terms within the 
brackets contain elements in common with the Poisson equation, for which sophisticated 
methods of solution exist (e.g., expansion by spherical harmonics and FFT convolution). 
When /(v) is a Maxwellian, a harmonic expansion gives 



F = -Al 2 (l + ^j) G(lv) , 



(21) 



and 



v 



${lv) - G(lv) 



G(lv) , 



(22) 



(23) 



where D\\ (D±) is the component of the diffusion tensor parallel (perpendicular) to the 
particle's velocity in the co-moving frame, $ is the error function, and G(x) = $(x) — 
x&(x)/2x 2 . Since dv/dvi — Vi/v, the diffusion tensor fills out as 



D. 



a/3 



V a V{3 



+ 2^/3 Da 



(24) 



Equations (21) thru (24) present a fairly complete non-relativistic kinetic theory. 
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The time td required for a particle to be deflected by 90 degrees is a good measure of 
how long it will take a non-thermal injection to equilibrate (Spitzer 1962). This time is given 
by the condition 

D ± t, = v 2 . (25) 

Thus, for particles whose root mean square velocity is that of the group, v = y^3kT/m, we 
have lv = 1.225, D± = 0.7UA (m/s) 2 /s, and 

{3kT\ 3/2 1 

t s = ttt^t. ,x seconds . 26 

A slowing-down viscous time scale can be similarly defined. One of our primary goals is 
to find a relativistically correct version of Equation (26) — the majority of problems (even 
apparently difficult ones, such as nonthermal radiation from galaxy clusters) can be resolved 
with a glance at this simple formula. 

Chandrasekhar (1957) approached the problem quite differently. Consider replacing 
the truly discrete many-body potential $ = J2i e i e 'i/\ r ~ r 'l with the smoothed version 
$(r) = J ep(r')/|r — r'\d 3 r'. A distribution evolving only due to such a smooth, mean field 
force is 'collisionless' — i.e., no particle ever leaves its volume of (single particle) phase space. 
But in moving to an integrable charge distribution, we have lost track of small fluctuations 
in the number of particles in a given neighborhood. 

Accumulated over a small time, these fluctuations cause a random kink in a particle's 
path. Alone, this kink diffuses the probability of finding the particle with a particular ve- 
locity until all velocities are equally likely. Since, however, a Maxwellian distribution for 
this probability inevitably sets in, Chandrasekhar reasoned that an associated viscous term 
must exist to bring diffusing particles back to the group, so that individual particles move 
through phase space according to the equations of motion 

£-* < 27 > 

^ = - + F d + Q.r(t), (28) 
at m 

where F is a smooth mean field term (plus any external forces) , Fd is dynamical friction that 
occurs when more particles hit the front end of a moving particle than the back, Q is the 
'square root' of the diffusion tensor, = QikQjk, and T(t) is a trivariate Gaussian random 
vector with 

(r 4 (f)> = , (29) 
(T t (t)T J (t)) = 5 tJ 5(t-t'). (30) 
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We show in Appendix A why Equations (2) and (28) are equivalent. 

Equations (27) and (28) reappropriate kinetic theory as a dynamics that is conceptually 
different than Hamilton's: Analytical dynamics follows the continuous motions of each of 
N particles in a 2 N- dimensional phase space; it conserves energy, maintains the number of 
particles in a volume of phase space as a constant, and is reversible. Analytical dynamics 
has been described as "a gradual unfolding of a contact transformation." 

By contrast, statistical dynamics follows the motion of particles accumulated over some 



coarse time At — indeed, since y |Av| 2 /At — > oo as At — > 0, statistical dynamics is not 
even defined continuously. Statistical dynamics ignores all but the 6 phase space coordi- 
nates of the particle in question, representing the effect of the other degrees of freedom as 
friction and noise — these terms cause a flux of particles through a surface of the single par- 
ticle phase space. Statistical dynamics does not conserve energy, having at its heart the 
dissipative force Fd, yet this force is precisely what is required to bring the assembly as a 
whole to a Maxwellian distribution characterized by a constant energy. Statistical dynamics 
is irreversible. Chandrasekhar described it as the "gradual unfolding of a Markoff chain." 

Solving Equation (1) by moving many particles according to the equations of motion 
(27, 28) is sometimes called a 'Monte Carlo' approach to kinetic theory. It is important to 
stress the historical development here to avoid the tincture of a computational trick that 
this name implies: the Langevin Equation (28) is equivalent to the Fokker-Planck equation, 
and implies no further approximation provided collisions are short (see Eq. 30). It can 
be developed on physical grounds entirely separate from the Boltzmann equation, and in 
fact was developed some 14 years before a successful solution of the Fokker-Planck equation 
appeared. 

While we will use Equations (27) and (28) to solve the evolution of the distribution, we 
return to Equation (2) to derive the relativistically correct kinetic coefficients. The crucial 
distinction between non-relativistic theory and the relativistic version that follows is that 
the function 



is not a relativistically invariant 4-scalar, since in different coordinate systems different 
particles will be considered to be located simultaneously in a given volume. The integrals 
N = J f(t, r, p) d 3 p, and i = J vf(t, r, p) d 3 p, however, do form a 4-vector: 





(31) 




(32) 



Here d 3 p/e is covariant, and if i k is covariant, than the phase space distribution / must also 
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be covariant. Rewriting Equation (3) in the form 

ee'wff'{d 3 p/e){d 3 p'/e') d 3 x , 



(33) 



we reason that, since the phase-space density, d 3 p/e, and i k are covariant, ee'w must also be 
invariant. It follows that the integral 



W 



ki 



1 



ee / q k q l v rd da , 



(34) 



corresponding to Equation (12), must form a symmetric 4-tensor. Here v ret is the velocity of 
one particle in the rest frame of the other. Following a derivation similar to (10-13), we find 



B af} = W al3 /ee' = A- 



v' 2 5 a > ~ ''I.'' 



,/3 



(35) 



to be the desired manifestly covariant generalization of normal kinetic theory. 



It is reasonable to neglect the zeroth component of (35). The reason for this is that the 
change in particle energy q° is of second order with respect to small scattering angle, and 
thus W 00 and W 0a are of third or fourth order. Yet the Fokker-Planck equation is accurate 
only as far as second-order quantities. Leaving only the components of Equation (8) corre- 
sponding to momentum flux, and re-expressing B a p in a convenient form, we arrive at 



C = 



_d_ 

8u r 



where 



D al! = A J Z(u,u')/(u')A', 



f=-aj: 



and where the kernel is now given by 



d 



du 



-Z(u,u') 



Z(u') dV , 



(36) 

(37) 
(38) 



Z(u,u') 



w 5 a(3 - u a Ufi - u' a u'p + r(u a u' p + u' a up) 



(39) 



Here, u is the ratio of the momentum to the rest mass m . Also, 7 = a/1 + u 2 /c 2 , 7' = 

r = 77' - u . u'/c 2 , (40) 

w = C \lr 2 - 1 . (41) 



a/1 + u'7c 2 , 
and 
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In the non-relativistic limit, r — > 1, w — > |u — u'| and u — > v. Thus relativistic kinetic 
theory, which begins as manifestly covariant, also transitions to the non-relativistic limit at 
low velocity. From now on, when we refer to 'velocity' in a relativistic context, we mean 
u = p/rriQ. 

Equation (36) is valid as long as bremsstrahlung is not so dominant as to make particle 
collisions inelastic, thus affecting the cross section. Provided this is the case, it is simple 
enough to calculate the bremsstrahlung losses as an integral over the electron distribution. 

Equations (37) and (38) may be solved with direct numerical integration. The Langevin 
Equation (28) can then be used to update a large number of phase-space points, rather than 
obtaining an implicit solution of the differential Equation (6). Besides the benefits of speed 
and the guaranteed conservation of particle number, this approach also allows various new 
acceleration mechanisms to be incorporated: as an example we include diffusion due to 
Alfvenic turbulence and radiative losses within a single-species electron plasma in §§3 and 4. 
It remains to show that the approach conserves energy and that the energy equipartitions 
correctly. 



3. Comparison with Existing Theory 

Let us compare this result to that of Nayakshin and Melia (1998, NM98). Since NM98 
construct a kinetic theory around the unitless energy E = 7 — 1, while covariant theory 
measures the friction and diffusion coefficients in the generalized velocity u = 7V, direct 
comparison is not obvious. For example, in NM98 the energy exchange and viscous coeffi- 
cients are equivalent, whereas covariant theory finds the energy exchange by integrating the 
flux of energy through a volume of phase space. Covariant theory gives diffusion coefficients 
both parallel and perpendicular to a particle's velocity; the diffusion coefficients of NM98, 
measured in the scalar 7, do not. 

While the two theories measure different quantities, we may nevertheless place them in 
the same system of units for direct comparison. The unitless diffusion coefficient is 

D aa = ^J Z{u, «')/(«') ^u l2 du' , (42) 

where the c 2 term gives the diffusion coefficient in units of 1/s, as in NM98 (the kernel Z is 
measured as 1/c, the constant A as c 3 /s). With a change of the variable of integration, 

D aa = — Z{u,u')f{i)—di, (43) 
c J 7 
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the function 

D(u,u') = ^Z(u,u')^ (44) 

is now the equivalent of Equation (35) in NM98. As we shall see below, the difference 
between these two coefficients appears to be the reason why we do not recover Blasi's (2000) 
results, since the method of NM98 is only valid in the highly relativistic regime (see below). 

At first, the most glaring discrepancy is the absolute scale of the two theories: covariant 
theory (along with nonrelativistic Rosenbluth potentials) is divergent (Fig. la), while that 
of NM98 (Fig. lb) is not. 

There are also physically grounded discrepancies. We show in the appendix that the 
Fokker-Planck and dynamical friction approaches are equivalent. In the latter theory, it 
is physically clear that all particles must be slowed by the viscous term, since a moving 
particle is always struck more frequently on its leading side, regardless of what relation its 
speed has to the rest of the group. To use an analogy with another 1/r potential: a star 
is never accelerated by gravitational drag. Yet the energy exchange term in NM98 (Fig. 
2) changes sign, drawing all particles to the average (scalar) speed. In this case, one must 
multiply the frequency of collisions by the fractional change in energy, which allows particles 
to both gain and lose energy as the result of multiple collisions. Thus, the covariant friction 
is positive-definite, while the 'friction' term 0(7,7') °f NM98 is not. 

Perhaps the best way to compare the two theories, and to gauge their domain of validity, 
is to observe their predictions for the temporal evolution of a distribution. Crucially to the 
application of Section 4, we find that the theory of NM98 may not be applied in the extreme 
non-relativistic limit. 

We follow the description of NM98, beginning with a Maxwellian distribution and solv- 
ing the temporal evolution of the distribution via the Chang-Cooper implicit method rec- 
ommended by Petrosian and Park (1996), for a variety of transrelativistic temperatures 
6 = kT/m e c 2 (so that 9 ~ 10 is roughly 5 x 10 10 K). A fully relativistic distribution {6 ~ 5) 
remains in equilibrium, conserving energy and particle number to around 1%. The kinetics 
of NM98 are therefore robust for highly energetic particles. 

As the temperature falls near the electron rest mass energy, however, the distribution 
deviates by between 3 and 15%. At these intermediate temperatures ~ 1, an algorithm 
which explicitly conserves particle number (as suggested in NM98 Appendix A) is required 
to resolve the coefficients. Neither the stochastic nor the implicit algorithms maintain a 
transrelativistic distribution in equilibrium. 

At 6 ~ 0.1, the distribution ceases to be relativistic, the function exp(— 7/6*) suffers from 
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underflow errors, and the kinetics of NM98 do not maintain an equilibrium distribution. In 
Blasi (2000), the distribution prior to heating lies at — 0.01; it is therefore the failure 
of NM98 to correctly transition to the nonrelativistic limit that leads to the inaccuracy in 
Blasi's (2000) formulation of the cluster's nonthermal emissivity. 

4. A Solution of the Equations in Covariant Kinetic Theory 

Figure 3 shows the kinetic coefficients calculated for a Maxwellian distribution 



at a temperature 10 8 K, compared to the known analytic versions given by Equations (21- 
23). It's important to stress here that, while most of our results are presented as functions 
of the magnitude of the velocity, this quantity plays no role in the calculation itself. 

To make this figure, a large number (here, one billion) of particles are first initialized, 
scaled, and interpolated onto a grid, generally using the cloud-in-cell or triangle-shaped- 
cloud interpolants well-known from many-body mesh calculations. This three dimensional 
grid represents a discrete version of the distribution /(v') on the right hand side of Equations 
(15) and (16). For any velocity v', a subfunction interpolates values off the grid, and the 
result is an apparently continuous function. 

We now concern ourselves with a single value of v — the left hand side of Equation (15)— 
and carry out a Monte Carlo integration of the kernel (Eq. 10) using the Vegas algorithm of 
the Cuba integration library. A maximum of 2500 function evaluations, beginning with 300 
subdivisions and adding 500 new subdivisions in a refinement step gives the best speed-to- 
accuracy ratio. For this one particular x, y, and z velocity, we sum over f3, resulting in six 
independent diffusion coefficients and three force coefficients comprising a structure. 

We do this for some number of v values, creating a second grid, each point of which is 
a structure containing the six D a p or three F a . Two computational points are relevant here: 
first, we perform the sum over (5 within the function call, rather than integrating three times 
and adding the answers. Second, the algorithm is about ten times faster when a vector of 
many functions is passed to the integrator in a single integration call, rather than defining 
the function and calling the integrator many separate times. A lookup table is created so 
that in calling the integrator to function number one, say, the code automatically knows that 
it's looking for the D xx value for v=(lx 10 8 , -2 x 10 8 , 2 x 10 8 ) m s -1 . 

Again, so long as the grid is sufficiently fine, a subfunction interpolates D a p and F a off 




(45) 
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the grid for any v. We now have an apparently continuous function solving (Eqs. 19-20). 

The f-grid points are divided evenly among several processors using MPI. Each process 
accumulates its local distribution, which is then summed and re-broadcast to all processes. 
A process solves the v values belonging to it, and finally the full grid is gathered. In this 
fashion, one can avoid the passing of particles. 

Finally, for display purposes, we find the magnitude of each particle's velocity and 
friction, and we diagonalize its diffusion tensor. That is, the coefficients in Figure 1 play no 
role in evolving the distribution; the solution makes no isotropic assumption and is naturally 
calculated in the lab frame. 

To evolve the distribution, each of the particles is updated using a first-order stochastic 
equation solver. Second order schemes exist (Qiang 2000), but the multiple function eval- 
uations are prohibitively expensive. The tensor Q is found by diagonalizing D via Jacobi 
rotation, taking the square root of the diagonal components, and then transforming back. 

Making a transition to a relativistic Maxwellian plasma, we instead use the distribution 
function 

'<"> = ^kTcK^lkT) ^i-im^lkT) , (46) 

where K2 is the modified Hankel function of order 2. Again, we are working with u = p/mo, 
the relativistic generalization of the velocity, and we use the term 'velocity' exclusively with 
this meaning. With a substitution of f(u) in Equation (37), the resulting relativistic kinetic 
coefficients may be defined as 

D R (u) = (f>(u)D NR (u) , (47) 

F R (u)=^(u)F NR (u). (48) 

For ease of use, we have expanded (f>(u) and ip(u) — functions we call the enhancement fac- 
tors, these being the multiplicative differences between the non-relativistic and relativistic 
expressions — as 

3 

0(«) = ^a i « i , (49) 


and similarly for if)(u), where again the coefficients aj are expanded as functions of temper- 
ature, 

°- = P>{t^k)'- («» 

The resulting coefficients are given in Table 1; the polynomial functions over bj fitted to 
each cij are shown in Figure 4. In Figure 5 we show the actual enhancement factors <f>(u) and 



14 



ip(u), together with the polynomials reconstructed from Table 1. Thus, with a few numbers, 
we can readily calculate the relativistic kinetic coefficients. 

These enhancements can now be used to redefine the equilibration and viscous time 
scales, one of our primary goals in this transition to a covariant kinetic theory. We have 



where the average velocity is evaluated from (7 — l)m e c 2 = (3/2)akT. The factor a is 1 
in the non- relativistic regime, but approaches 2 when 7 >> 1. We may estimate it as the 
piecewise function 



Now the length of time required for a relativistic distribution to relax to equilibrium may be 
measured against that for a non-relativistic distribution by the ratio t^n/t^NR = v 2 /u 2 (f>(u), 
which is plotted in Figure 6. (Again, the definition of the Spitzer time is valid at the average 
scalar velocity of the group. The Spitzer time is not particularly useful for calculations: 
rather, it is designed to give an intuitive grasp of the amount of time needed for equilibration, 
and we give its relativistic generalization as motivation for the re-examination of a wide class 
of problems.) Note that a fully relativistic plasma equilibrates in l/100th of a Spitzer time, 
t s - 

The evolution of the non-relativistic and relativistic distributions — at 10 8 and 2.5 x 10 9 
K, respectively — is shown in Figures 8 and 9, as a function of the Spitzer time, t s . The 
non-relativistic distribution requires about two Spitzer times to equilibrate, whereas the 
relativistic one reaches equilibrium in only 0.2 t s . 

When the energy is constant, or when the injection rate is known, we use normalization 
coefficients to maintain constant energy in the distribution; an example of how these coef- 
ficients vary with time is plotted in Figure 10. Energy is conserved to an accuracy of 10~ 2 , 
about the same as the precision of our integration. 

However, when energy is not constant, this simple procedure is not feasible, and we 
must increase the grid size to a minimum of 64 3 and the maximum number of function 
evaluations to 10,000. In this case, normalization is no longer required. For example, in 
Figure 11, we allowed the distribution to relax for a full Spitzer time before turning on 
stochastic acceleration, conserving energy all the while. 



D± R t diR = u 



(51) 




(52) 
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5. Sample Application to Cosmic Ray Acceleration in Galaxy Clusters 

Evidence for the presence of relativistic electrons in the intracluster medium (ICM) 
is provided primarily by diffuse synchrotron emission at radio wavelengths. Given the ~ 
Mpc-size of the synchrotron features, and the fact that the radiative lifetime of the electrons 
appears to be orders of magnitude shorter than the time required to cover such distances, 
the presence of this emission also suggests that electrons are being accelerated out of the 
ICM thermal pool. 

Typically, it is assumed that the radiating particles must be continuously re-accelerated 
on their way out — these are 'primary' models (Jaffe 1977; Schlickeiser et al. 1987, Brunetti 
et al. 2001, Ohno et al. 2002). However, protons could diffuse throughout the cluster 
without radiating, all the while colliding with thermal ICM ions and cascading into the 
observed relativistic electrons — these are 'secondary' models (Dennison 1980). Secondary 
models relieve the difficulties of continuously accelerating electrons against radiative and 
Coulomb cooling, however Blasi and Colafrancesco (1999) suggest that protons could not 
provide enough secondary electrons to describe the radio emission without also producing 
7T° gamma decays in excess of the number observed by EGRET. 

Coma's radio emission is accompanied by hard X-rays of similar spatial distribution; 
this and its regular morphology classify the Coma as a 'radio halo'. In contrast, 'radio relics' 
are typically irregular and concentrated toward the cluster's periphery (e.g., Feretti 2003). 

Interpreted as thermal bremsstrahlung, the observed X-rays imply an ICM tempera- 
ture between 8 and 9 keV. Coma's X-ray emission cannot be described entirely as thermal 
bremsstrahlung, however. The Rossi X-Ray Timing Explorer (Rephaeli and Gruber 2002) 
and BeppoSAX (Fusco-Femiano et al. 2004, hereafter FF04) have each made two observa- 
tions of the Coma, both claiming to find a hard X-ray tail (HXR) in excess of a thermal flux. 
This is a controversial claim: the initial analysis of Fusco-Femiano et al. (1999, FF99) was 
flawed — one of three spectra was counted twice — leading to a 2.5 a overstatement of the 
confidence level, according to Rossetti and Molendi (2004, RM04). RM04 flatly state that 
the second BeppoSAX observation shows no evidence of a hard tail. FF04, however, find a 
significant drop in the x 2 value for the fit to a hybrid thermal-power-law model (x 2 = 1.20 
for 7 dof) as opposed to the purely thermal model (x 2 = 4.10 for 9 dof). This improvement 
cannot be matched by a two-temperature model, since the second temperature (kT > 50 
keV) is considered unrealistic. 

Rephaeli (1979) had predicted such a tail must exist when the synchrotron-emitting 
electrons inverse Compton scatter with the cosmic background. In principle, this origin could 
be used as a probe for the ICM magnetic field. Yet, assuming that both the synchrotron 
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radio and Compton HXR are produced by the same population of electrons, FF04 infer 
a magnetic field of B ~ 0.2 /xG — more than a factor of ten below the results inferred by 
Faraday rotation (B ~ 6.0 /xG). 

For this reason, nonthermal bremsstrahlung (Blasi 2000) was proposed as an alternative 
emission mechanism to inverse Compton. Here, nonthermal electrons are directly accel- 
erated and directly radiate via e — e bremsstrahlung. Petrosian (2001) has argued that 
bremsstrahlung is not sufficiently efficient to cool the distribution before an unacceptably 
large amount of energy is given to the ICM, shifting the thermal body of bremsstrahlung 
above its observed flux. Because this problem is simple (it requires only stochastic acceler- 
ation, bremsstrahlung and Coulomb cooling), and because it straddles the transrelativistic 
regime which our theory handles uniquely, we have chosen it as an example of the technique's 
power and simplicity. 

Our calculations agree with the analysis of Petrosian. They suggest that the source of 
confusion is the inability of NM98 and DL89's kinetics to reduce to the nonrelativistic limit. 
We find that the stochastically gained energy heats the body of the distribution, not just 
the tail, on a Spitzer timescale of some tenths of a Myr; compare this to the tenths of Gyr 
required for Blasi's model. Similarly, the Spitzer time requires that any merger event must 
have occurred within a few Myrs, a vanishingly small window of time. 

With this motivation in mind, let us state the problem more precisely. We use the 
pitch-angle averaged diffusion coefficient (Dermer, Miller, and Li 1996), 

c%p 2 AVA (r B k y- 2 P yp , (53) 

to represent the resonant interaction of particles with Alfven waves. Here p = 7/5, r B = 
m e c 2 /eB, t]a = 0.07 is the fraction of magnetic energy in Alfven waves, q = 5/3 is the 
Kolmogorov constant, L c = 10 kpc is the size of the galaxy, Pa = va/c (where va = 
B/ yj A.im p m p is the Alfven velocity), and k Q = 2n/L c is the largest-scale wavenumber. Our 
point of departure is that of Blasi's calculation (n = 4 x 10~ 4 cm -3 , B = 0.8 /xG, and T = 7.5 
keV), but we also look at a small range around these figures. 

In Figure 10, we show the evolution of the distribution, begun at 7.5 keV, with a window 
of 0.1 Spitzer times between each curve. One full Spitzer time elapses before stochastic 
acceleration from a 0.8 /xG field is turned on (corresponding to the first curve on the left). 
We see that rather than stochastic acceleration producing a high-energy tail, the high rate 
of thermalization provides energy to the whole distribution. 

The emissivity of this distribution is calculated self-consistently using particle number 
distributions and electron-proton bremsstrahlung cross section, 



D(u) 



1 



q(q + 2) 



-17- 



dE 



dV dt dv 
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where da/dk is the relativistic bremsstrahlung cross section averaged over all angles, and 
expanded to the (p ) 6 th order. Including the Elwert correction factor, this is (Haug 1997) 
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where a = e 2 /hc is the fine structure constant, r Q = e 2 /m e c 2 is the electron radius, k = 
hu/m e c 2 is the photon energy in units of the electron mass energy, p Q = joPo = u Q /c is the 
initial electron momentum in units of m e c, e Q = 7„ is the total energy in rest mass units, 
a = ae /p , and 



Pf= [k 2 -p -2k^l-p^j 



1/2 



(56) 



In the energy regime considered, the error of this expansion is a few percent — the contribution 
of e~-e _ and e~-e + bremsstrahlung is of this order, so we exclude them from consideration. 

Once the cross section in Equation (55) is found, we calculate the observed X-ray flux F x 
under the assumption of constant density and use the accepted volume V and distance dh for 
Coma to write F x = V j / \47id 2 L ) , assuming a (dimensionless) Hubble constant h = 0.6. The 
calculated spectrum is shown in Figure 11, together with the BeppoSAX data, at spacings 
of 0.8f 8 . 

For completeness, we also show the evolution for nine combinations of electron density 
and magnetic field. In the bottom left corner — higher density and lower magnetic field — a 
nonequilibrium tail never heats. In the top right — lower density and higher magnetic field— 
the soft X-ray flux immediately moves above the observed limits. In none of these figures 
do we address the other free parameter, the initial temperature. And, while the calculation 
is done via the relativistically correct kernel, in no case is a relativistic tail produced before 
the nonrelativistic body heats to some excessively high energy. 

All of these figures evolve over At s = 2 Myr, some factor > 10 2 below that reported by 
Blasi. The corrected Spitzer timescale is irreconcilable with the HXR observations (Fig. 11), 
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as in no case is a nonthermal tail accelerated before the ICM equilibrates the added energy 
and reaches an unacceptably high temperature. This is consistent with the result reported 
by Petrosian (2001): the erstwhile lack of a kinetic theory that is applicable throughout the 
transrelativistic regime was the source of discrepancy. 

6. Concluding Remarks 

We have described a novel method for solving the covariant Landau equation, uncon- 
strained by assumptions of isotropy or low particle energy. Using simple polynomial fitting 
formulae for the kinetic coefficients, we have described an efficient numerical technique for 
determining the temporal evolution of an arbitrarily relativistic particle distribution, which 
may also be subject to energy injection from the anomalous acceleration of particles, e.g., in 
shocks or scattering with Alfven waves. We anticipate that this method will find widespread 
application not only in astrophysics, but other physics disciplines as well. 

To demonstrate the power of our technique in understanding the behavior of high- 
energy plasmas in astrophysics, we have examined one of the currently considered models — 
bremsstrahlung emission by a nonthermal tail — for the hard X-ray emission in galaxy clus- 
ters, such as Coma. Our conclusion is that the time required by the underlying plasma to 
attain equilibrium is far too short compared to the stochastic acceleration time scale for 
any nonthermal high-energy extension to survive longer than ~ 0.5 Myr. It seems unlikely, 
therefore, that cluster mergers could be responsible for energizing the plasma turbulence 
required to sustain the nonthermal particle emissivity. 

The application of covariant theory to relativistic electrons in the ICM remains highly 
relevant. As we suggest here, the time required for relaxation may vary by orders of magni- 
tude from the currently used quantity. In a future paper, we will examine the production of 
7T° gamma rays and the associated cascade of charged leptons during proton-proton collisions 
in the ICM. In addition, we have not yet completely resolved the complete set of differences 
between covariant theory and the kinetics of DL89 and NM98, which amounts to better 
defining the range of validity for the latter approximate treatments. We are performing a 
detailed comparison of their energy exchange and temporal evolution properties, and these 
results too will appear in a future paper. 
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8. Appendix A: Equivalence of the Langevin and Fokker Planck approaches 

We argue following Zwanzig (2001). Consider the stochastic equation 



dv 
~dt 



F d (v)+A(t) , 



(57) 



where the noise A(t) is Gaussian, with zero mean and delta-correlated second moment, 

(A(t)A(t)) = 2B8(t - t') . (58) 



Let us now find the probability distribution f(v,t), averaged over the noise. Now, f(v,t) is 
conserved: 

f) f r) I rlii \ 

. (59) 



dt dv \ dt 

Thus, substituting dv/dt we arrive at the stochastic differential equation 



dl = _ d_ 

dt dv 



F d (v)f(v,t)+A(t)f(v,t) 



-Lf-l v Mt)f, 



(60) 



where L$ = (d/dv)F d (v)§ is an appropriately defined operator. The solution to this equa- 
tion is 

dse-^ L —A(s)f(v,s). (61) 
Placing Equation (61) back into Equation (60), we arrive at 

= -Lf(v, t) - ^Mt)f(v, 0) + ^A(t) ^ dse-^ L ^A(t)f(v, s) , (62) 



dt 



which is readily observed to be the Fokker Planck equation, 

!</(«.«)> =-I^(«)(/(«.«)>+|b|(/(«.«)> 



(63) 



once the average of f(v , t) is taken, and the zero average and unit correlation of A is imposed. 
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Table 1: Expansion coefficients (for Eqs. 46 and 47). 
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Fig. 1. — A comparison of our theory with that of Nayakshin and Melia 1998. Our kernel 
(Eq. 44) is divergent when 7 = 7', while that of NM98 is not. Taking 7' = 1.05 to be 
effectively nonrelativistic, we can see that the theory of NM98 does not duplicate the results 
of Rosenbluth (1957). 
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Fig. 2. — NM98's viscous coefficient 0(7,7') is n °t positive-definite, as it is in standard 
kinetic theory (Fig. 3). 
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Fig. 3. — The theory of NM98 breaks down in the non-relativistic limit, perhaps due to 
limited machine precision as the Lorentz factor approaches unity. Here, a distribution begun 
at equilibrium (dashed) evolves to 1/4 (dash-dotted) and 4 (solid) Spitzer times under the 
influence of NM98's kinetic coefficients. The quantity 9 is the fraction of thermal energy to 
electron rest mass energy; 9 = 10 is roughly 5 x 10 10 K. 




Fig. 4. — Analytic (plus signs) and calculated (circles) D(v) and F(v). The Viscous Force 
and the Diffusion coefficients (perpendicular on top, parallel below) are normalized to one. 
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Fig. 5. — The four parameters e^, and the fitting polynomials (shown as smooth continuous 
curves) bj (T / 5 x 10 8 K) J . We have expanded the kinetic coefficients for a particular 
temperature — and, because this expansion changes as the temperature is raised, we again 
expand each of the coefficients as though they were polynomial functions of T. We present 
this figure not simply to confirm the accuracy of the expansion, but so that a rough estimate 
can be made by reading off the a; for the temperature of interest. 
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Fig. 6. — Enhancement factors <f>(u) and ip(u), reconstructed from Table 1, together with 
the raw data. 




Fig. 7. — Ratio of relativistic to non-relativistic relaxation time scales, as a function of 
temperature. 
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Fig. 8. — A distribution of particles (dotted) initialized with equally likely v x ,v y , and v z , 
whose standard deviations are that of a Maxwellian with temperature T = 10 8 K, relaxes 
(to the analytic form, shown as a solid curve) in a time 2.1 t s . The top four frames are in 
multiples of 0.1t s , the remaining 8 of 0.2t s 
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Fig. 9. — An initially cold distribution of particles (dotted) with the same energy as a 
relativistic Maxwellian at T = 2.5 x 10 9 K (thick solid) relaxes in 0.2 t s . Top four frames are 
in multiples of 0.015i s , remaining 8 of 0.025t s . For comparison, we show the non-relativistic 
distribution (Eq. 42; thin solid) at this temperature. 
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Fig. 10. — For a 32 3 grid, the kinetic coefficients are correct only to 2 significant figures; since 
they are self-consistent even this small deviation can lead to run-away heating, and we must 
"normalize" with a slowly evolving coefficient, plotted here as a function of time in units of 
the Spitzer time t s . This slight deviation from strict energy conservation can be removed by 
using instead a 64 3 grid (with its accompanying 4 significant figures of accuracy) . 
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Fig. 11. — The distribution evolves self-consistently from an initial temperature of 7.5 keV 
(curve on the left), as it is heated by Alfvenic diffusion and cooled by Coulomb diffusion and 
bremsstrahlung emission (curves moving progressively to the right) as the entire body heats 
over a period of only 4t s = 2 Myr. 
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Fig. 12. — X-ray spectra produced by the distributions shown in Figure 10 are in panel e. 
The others are from a range of possible values of magnetic field and density. Each line is 
plotted at an integer of 0At s . In a,d,g, and e, we have allowed 8t s to pass, to show that 
the body of the distribution heats by the time a nonthermal tail is produced; all other 
figures evolve over 4t s . Data points (error bars not shown) are from the first observation of 
BeppoSAX (Fusco-Femiano et al. 1999). Dashed line is equilibrium flux at a temperature 
of 8.21 keV. Either the distribution is never driven from equilibrium (d,g,h, and i), the body 
receives an excessive amount of heat (e,f ) or the whole distribution goes over to a power law 
(a,b and c). 



